#!/usr/bin/env python

# -*- coding: utf-8 -*-

"""
	Copyright (C) 2008  Michael Ughetto <telimektar1er@gmail.com>

	This file is part of Danuve.
	Danuve is a easy to use numerical tool meant to solve chain of oscillators

	Danuve is free software: you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation, either version 3 of the License, or
	(at your option) any later version.

	Danuve is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with Danuve.  If not, see <http://www.gnu.org/licenses/>.
"""  

import sys, os
import threading

class oscillators_chain(threading.Thread):
	def __init__(self, masses, positions, speeds, \
			 ressorts_left, ressorts_right, delta_t, time):
		threading.Thread.__init__(self)
		print "Initializing the oscillator chain...\n"
		self.masses = masses
		self.positions = positions
		self.positions.append(0) # This is a virtual position
		self.speeds = speeds
		self.ressorts_left = ressorts_left
		self.ressorts_right = ressorts_right
		self.delta_t = delta_t
		self.time = time
		self.work_time = False

	def set_solve_info(self, output_file, algo='euler'):
		"""set_solve_info define non-physics parameters for object

		output_file -- file URI where to write the results
		algo -- the algo used for solving the system 

		"""
		self.output_file = output_file
		self.algo = algo
	
	def energy(self):
		pass

	def movement_equation(self, positions, i):
		"""movement_equation is used to determine the infinitesimal acceleration
		for the next step

		positions -- the list of positions
		i -- the position of the evaluated masse

		"""
		acc_1 = -1 * self.ressorts_left[i] * (positions[i] - positions[i - 1])
		acc_2 = -1 * self.ressorts_right[i] * (positions[i] - positions[i + 1])
		ds = (acc_1 + acc_2) / self.masses[i]
		return ds
		

	def euler_step(self):
		"""euler_step extrapolate the position of all masses after a little 
		delta_t by using the Euler Method.
		
		output_file -- file URI where to write the results
		
		"""
		future_positions = [0] * len(self.masses)
		future_speeds = [0] * len(self.masses)

		i = 0
		while i < len(self.masses):
			future_positions[i] = self.positions[i] + self.speeds[i] * \
				self.delta_t
			ds = self.movement_equation(self.positions,i)
			future_speeds[i] = self.speeds[i] + ds * self.delta_t
			self.output_file.write(str(future_positions[i]) + ' ')
			i += 1
		i = 0
		while i < len(self.masses):
			self.positions[i] = future_positions[i]
			self.speeds[i] = future_speeds[i]
			i += 1
		try:
			self.output_file.write('\n')
		except ValueError:
			pass

	def rk4_step(self):
		"""rk4_step extrapolate the position of all masses after a little 
		delta_t by using the Runge Kutta 4 method.
		
		output_file -- file URI where to write the results

		"""
		# No need midfuture speeds, for they don't count for the movement
		# equation.
		midfuture_positions_1 = [0] * (len(self.masses) + 1)
		midfuture_positions_2 = [0] * (len(self.masses) + 1)
		nearfuture_positions = [0] * (len(self.masses) + 1)

		k1x = [0] * len(self.masses)# I create all the k now,
		k2x = [0] * len(self.masses)# it make the calcul code clearer
		k3x = [0] * len(self.masses)
		k4x = [0] * len(self.masses)

		k1s = [0] * len(self.masses)
		k2s = [0] * len(self.masses)
		k3s = [0] * len(self.masses)
		k4s = [0] * len(self.masses)

		i = 0
		while i < len(self.masses):
			k1x[i] = self.speeds[i] * self.delta_t
			k1s[i] = self.movement_equation(self.positions, i) * self.delta_t
			i += 1
		i = 0
		while i < len(self.masses):
			k2x[i] = (self.speeds[i] + k1s[i] / 2.) * self.delta_t
			i += 1
		i = 0
		while i < len(self.masses):
			midfuture_positions_1[i] = self.positions[i] + k2x[i] / 2.
			i += 1
		i = 0
		while i < len(self.masses):
			k2s[i] = self.movement_equation(midfuture_positions_1, i) * \
				self.delta_t
			i += 1
		i = 0
		while i < len(self.masses):
			k3x[i] = (self.speeds[i] + k2s[i] / 2.) * self.delta_t
			i += 1
		i = 0
		while i < len(self.masses):
			midfuture_positions_2[i] = self.positions[i] + k3x[i] / 2.
			i += 1
		i = 0
		while i < len(self.masses):
			k3s[i] = self.movement_equation(midfuture_positions_2, i) * \
				self.delta_t
			i += 1
		i = 0
		while i < len(self.masses):
			k4x[i] = (self.speeds[i] + k3s[i] * self.speeds[i]) * self.delta_t
			i += 1
		i = 0
		while i < len(self.masses):
			k4x[i] = (self.speeds[i] + k3s[i] * self.speeds[i]) * self.delta_t
			i += 1
		i = 0
		while i < len(self.masses):
			nearfuture_positions[i] = self.positions[i] + k4x[i]
			i += 1
		i = 0
		while i < len(self.masses):
			k4s[i] = self.movement_equation(nearfuture_positions, i) * \
				self.delta_t
			i += 1

		# Now we have all k, we can go one !
		future_positions = [0] * (len(self.masses) + 1)
		future_speeds = [0] * (len(self.masses) + 1)

		i = 0
		while i < len(self.masses):
			future_positions[i] = self.positions[i] + (1 / 6.) * \
				(k1x[i] + 2 * k2x[i] + 2 * k3x[i] + k4x[i]) 
			future_speeds[i] = self.speeds[i] + (1 / 6.) * \
				(k1s[i] + 2 * k2s[i] + 2 * k3s[i] + k4s[i]) 
			self.output_file.write(str(future_positions[i]) + ' ')
			i += 1
		i = 0
		while i < len(self.masses):
			self.positions[i] = future_positions[i]
			self.speeds[i] = future_speeds[i]
			i += 1
		try:
			self.output_file.write('\n')
		except ValueError:
			pass

	def rk4_var_step(self):
		former_step = self.step
		

	def run(self) :
		"""run will launch the one_step method till the end of the 
		specified integration interval (starting from 0). 
		
		output_file -- file URI where to write the results
		algo -- The algorithm used to solve the system
		
		"""
		self.Terminated = False
		if self.work_time:
			self.timer = threading.Timer(self.work_time, self.stop)
			self.timer.start()
		print "Solving the described system ..." 
		t = 0.
		if self.algo == 'euler':
			print "Using Euler algorithm ..."
			while t < self.time and not self.Terminated:
				self.output_file.write(str(t) + ' ')
				self.euler_step()
				t += self.delta_t
		elif self.algo == 'heun':
			print "Using Heun algorithm ..."
			while t < self.time and not self.Terminated:
				self.output_file.write(str(t) + ' ')
				self.rk4_step()
				t += self.delta_t
		elif self.algo == 'rk4':
			print "Using Runge-Kutta 4 algorithm ..."
			while t < self.time and not self.Terminated:
				self.output_file.write(str(t) + ' ')
				self.rk4_step()
				t += self.delta_t
		elif self.algo == 'var_rk4':
			print "Using a variable step Runge-Kutta 4 algorithm ..."
			while t < self.time and not self.Terminated:
				self.output_file.write(str(t) + ' ')
				self.rk4_step()
				t += self.delta_t
		if self.Terminated == False:
			self.stop()

	def stop(self):
		"""stop will end the thread properly.
		
		"""
		self.Terminated = True
		if self.work_time:
			self.timer.cancel()
		print "\nResults can be found in", self.output_file.name
